Create a Base Map

A male Kentucky Arrow Darter collected from Wild Dog Creek.

First, lets create a kentucky state basemap

state <- map_data("state")
county <- map_data("county")
#apsu_point <- data.frame("x" = -87.353069, "y" = 36.533654)
ky <- county %>% 
  filter(region=="kentucky")

ggplot() + geom_polygon(data = state, aes(x=long, y = lat, group = group),
                        fill = "white", color="black") + 
           geom_polygon(data = ky, aes(x=long, y = lat, group = group),
                        fill = "gray", color="black") + 
  coord_fixed(xlim = c(-90, -82),  ylim = c(36.5, 39), ratio = 1.2) + 
  xlab("Longitude") + ylab("Latitude") + ggtitle("Kentucky")

Amazing. Now we have the state of Kentucky, lets add some more data.

kad_co_streams<- sf::read_sf(("~/Desktop/aculley_KAD_master_folder/KAD_base_maps/KAD_base_maps/kadstreamdata/kad_stream_data.shp"))
kad_watershed<- sf::read_sf(("~/Desktop/aculley_KAD_master_folder/KAD_base_maps/KAD_base_maps/HUC8_CONUS/HUC8_US.shp"))
kad_local <- read.csv("/Users/alexisculley/Desktop/aculley_KAD_master_folder/KAD_base_maps/KAD_base_maps/kad_localities.csv")
sur_mine <-  sf::read_sf(("/Users/alexisculley/Desktop/aculley_KAD_master_folder/KAD_base_maps/KAD_base_maps/Shapefile/5072/2019_activeMining_5072.shp"))

we just added stream, watershed, KAD locality, and surface mining data.

Now lets practice getting data straight from the internet! Lets get a kml file from github that I added.

add .kml file from the web:

local.kml <- readOGR("https://raw.githubusercontent.com/avculley/more_maps/main/kadlocal.kml")
## OGR data source with driver: KML 
## Source: "https://raw.githubusercontent.com/avculley/more_maps/main/kadlocal.kml", layer: "KADlocalities.csv Events"
## with 12 features
## It has 2 fields

Lets make this file usable

local_points <- cbind(local.kml@data, local.kml@coords)
local_points[2] <- NULL
local_points[4] <- NULL
colnames(local_points) <- c("name","x","y")

x = ggplot()+ 
  geom_polygon(data = state, aes(x=long, y = lat, group = group), fill = "white", color="black") + 
  geom_polygon(data = ky, aes(x=long, y = lat, group = group), fill = "gray", color="black") +
    geom_point(data=local_points, aes(x = x, y = y, color = name), size = 4, alpha = 0.8) +
  coord_fixed(xlim = c(-90, -82),  ylim = c(36.2, 39.2), ratio = 1.2)

x

yay they worked!

Now we need to graph these to see what they look like. We added a datset including watershed outlines, localities, and streams.

ggplot(kad_co_streams)+geom_sf()

kad_co_rivers <- subset(kad_co_streams, NAME == "South Fork Kentucky River" | NAME == "Middle Fork Kentucky River" | NAME == "North Fork Kentucky River" | NAME == "Kentucky River")

KAD streams looks great! But.. Hold up, we cant just plot this because it will take 7 years since this dataset contains the entire United States. Lets filter the watershed data to only include the areas of Kentucky we care about. You can do this is a lot of ways, but I am just going to subset the dataset by ID’s I know from the dataset.

HUC8_KAD <- subset(kad_watershed, HUC8 == "05100201" | HUC8 == "05100202" | HUC8 == "05100203" | HUC8 == "05100204")
ggplot(HUC8_KAD)+geom_sf()

Now that we have data on the scale we need, we have proof they are actually into R.

Let’s plot them on top of each other. First, we need to set the crs for the datasets so they match.

st_crs(kad_co_streams)
## Coordinate Reference System: NA
kad_co_streams_crs = st_set_crs(kad_co_streams, 4326)
st_crs(kad_co_rivers)
## Coordinate Reference System: NA
kad_co_rivers_crs = st_set_crs(kad_co_rivers, 4326)
st_crs(sur_mine)
## Coordinate Reference System:
##   User input: NAD83(NSRS2007) / Conus Albers 
##   wkt:
## PROJCRS["NAD83(NSRS2007) / Conus Albers",
##     BASEGEOGCRS["NAD83(NSRS2007)",
##         DATUM["NAD83 (National Spatial Reference System 2007)",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4759]],
##     CONVERSION["Conus Albers",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",23,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-96,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",29.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",45.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Data analysis and small scale data presentation for contiguous lower 48 states."],
##         AREA["United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming."],
##         BBOX[24.41,-124.79,49.38,-66.91]],
##     ID["EPSG",5072]]
sur_mine_crs = st_set_crs(sur_mine, 4326)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
ggplot() + 
  geom_sf(data = HUC8_KAD, color = "black") +
  geom_sf(data = kad_co_streams_crs, color = "blue") +
  ggtitle("Kentucky Arrow Darter Historical Range") + 
  coord_sf()

Now that we can plot them on top of each other, lets try and make it look less awful. Lets clip the stream data by the watershed polygons.

kad_watershed_streams = st_intersection(kad_co_streams_crs, st_make_valid(HUC8_KAD))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
kad_rivers_streams = st_intersection(kad_co_rivers_crs, st_make_valid(HUC8_KAD))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
test1 = ggplot() + 
  geom_sf(data = HUC8_KAD, color = "black") +
  geom_sf(data = kad_watershed_streams, color = "blue") +
  ggtitle("Kentucky Arrow Darter Historical Range") + 
  coord_sf()

test1

Wow! Now we have all the streams in the watersheds the Kentucky Arrow Darter are found in. This doesn’t tell us a lot though does it?

localmap = ggplot() + 
  geom_sf(data = HUC8_KAD, color = "#CC79A7", fill = "white") +
  geom_sf(data = kad_watershed_streams, color = "#56B4E9") +
  geom_sf(data = kad_rivers_streams, color = "#0072B2") +
      geom_point(data=local_points, aes(x = x, y = y, color = name), color="black", size = 4, alpha = .8) +
      geom_text(data=kad_local,aes(x=long,y=lat,label=name), color="white", vjust=1.6, hjust=-.11, size=2.5, fontface="bold") +
  geom_text(data=kad_local,aes(x=long,y=lat,label=name), color="black", vjust=1.5, hjust=-.1, size=2.5, fontface="bold") +
  labs(x="Longitude",y="Latitude", title="Kentucky Arrow Darter Localities", fill = "No. of Stations") +
  north(kad_watershed_streams, location = "bottomleft", scale = 0.05, symbol = 12) +
  theme(legend.position = "bottom")+
    annotation_scale(location = "bl", width_hint = 0.5) +
  theme(legend.position = "right") +
  coord_sf()
  
#other north arrows
#annotation_north_arrow(location = "br", which_north = "true", 
       # pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
       # style = north_arrow_fancy_orienteering)
#other kad locality file
#geom_point(data = kad_local, aes(x = long, y = lat), color="black", size = 3, alpha = 0.8) +

localmap

Now time to make the inset map.

state <- map_data("state")
county <- map_data("county")
#apsu_point <- data.frame("x" = -87.353069, "y" = 36.533654)
ky <- county %>% 
  filter(region=="kentucky")

ken_ref_map = ggplot() + 
  geom_polygon(data = state, aes(x=long, y = lat, group = group),
                        fill = "white", color="black") + 
           geom_polygon(data = ky, aes(x=long, y = lat, group = group),
                        fill = "gray", color="black") + 
  geom_sf(data = HUC8_KAD, color = "#CC79A7") +
  coord_sf(xlim = c(-90, -82),  ylim = c(36.5, 39)) + 
  theme(panel.background = element_rect(fill = "lightblue"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
        axis.line=element_blank(), axis.text.x=element_blank(), axis.text.y=element_blank(),axis.ticks=element_blank(), 
        axis.title.x=element_blank(), axis.title.y=element_blank())
ken_ref_map

Now lets combine our kad locality map with our inset map!

gg_inset_map1 = ggdraw() +
  draw_plot(localmap) +
  draw_plot(ken_ref_map, x = 0.05, y = 0.65, width = 0.3, height = 0.3)

gg_inset_map1

leaflet(kad_local) %>% 
  addTiles()%>%  
  addCircleMarkers(data = kad_local,
             lat = ~lat, lng = ~long, label = kad_local$name,
                   weight = 2,
                   color = "grey",
                   fillColor = kad_local$name,
                   fillOpacity = 0.7,
             stroke = FALSE, opacity = 0.3)%>%
  addWMSTiles("https://basemap.nationalmap.gov/arcgis/services/USGSTopo/MapServer/WmsServer", layers = "0") %>%
  addMiniMap(zoomLevelOffset = -4) %>%
  addScaleBar()

Bog Turtle Home Range Example

For this example you are going to obtain records for Glyptemys muhlenbergii) (Bog Turtle) within the eastern United States.

Description: The bog turtle is the smallest turtle in North America, rarely exceeding three or four inches in length and weighing only about four ounces. Its orange to yellow patch on either side of the neck easily distinguishes it from other turtles. Bog turtles emerge from their muddy hibernation in early to mid-April and by early May are actively seeking a mate. Adults are sexually mature at five to eight years of age. In June or July, the female lays a clutch of one to six small white elliptical eggs in a shallow “nest” she digs in a clump of sphagnum moss or tuft of grass above the water line. After seven or eight weeks of being incubated by the sun, the inch-long hatchlings emerge. Because they are born so late in the year, the hatchlings often spend their first winter near the nest.

Habitat: Bog turtles live in the mud, grass and sphagnum moss of bogs, swamps, and marshy meadows. These wetlands are usually fed by cool springs flowing slowly over the land, creating the wet, muddy soil needed by the turtles.

Range: There are two distinct populations of the bog turtle separated by about 250 miles. The northern population is found from New York

The current distribution of the Bog Turtle according to the FWS.

turt <- gbif("glyptemys", species = "muhlenbergii", ext = c(-130,-70,20,60),
                   geo = TRUE, sp = TRUE, download = TRUE,
                   removeZeros = TRUE)

The Bog turtle has 146 observations across the eastern coast, the next step is to prepare it for display. To do this, we are going to create an x,y dataset, obtain a base map of the US, and create our map of known Bog turtle occurrences.

turt_df <- cbind.data.frame(turt@coords[,1],
                                  turt@coords[,2])

colnames(turt_df) <- c("x","y")

us <- map_data("state")

ggplot(data = turt_df, aes(x=x, y=y)) +
  geom_polygon(data = us, aes(x=long, y = lat, group = group),
               fill = "white", color="black") +
  geom_point() + xlab("Longitude") + ylab("Latitude") +
  coord_fixed(xlim = c(-84,-72), ylim = c(34.8,46)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Bog Turtles in the Eastern US") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "lightblue"))

Let’s compare to the actual FWS data.

Honestly… not that bad.

Yay! Now let’s import some climate variables to create a model- I’m stealing the ones Dr. Gentry used in class.

bioclim <- getData(name = "worldclim", res = 2.5, var = "bio", path = "./")

names(bioclim) <- c("Ann Mean Temp","Mean Diurnal Range","Isothermality","Temperature Seasonality",
                    "Max Temp Warmest Mo","Min Temp Coldest Mo","Ann Temp Range","Mean Temp Wettest Qtr",
                    "Mean Temp Driest Qtr","Mean Temp Warmest Qtr","Mean Temp Coldest Qtr","Annual Precip",
                    "Precip Wettest Mo","Precip Driest Mo","Precip Seasonality","Precip Wettest Qtr",
                    "Precip Driest Qtr","Precip Warmest Qtr","Precip Coldest Qtr")

bio_extent <- extent(x = c(
  min(turt_df$x),
  max(turt_df$x),
  min(turt_df$y),
  max(turt_df$y)))

bioclim_extent <- crop(x = bioclim, y = bio_extent)
bioclim_model <- bioclim(x = bioclim_extent, p = cbind(turt_df$x, turt_df$y))
presence_model <- dismo::predict(object = bioclim_model, 
                                 x = bioclim_extent, 
                                 ext = bio_extent)

In this script bioclimatic variable were generated and descriptive names were provided. Because the bioclimatic variables are global, an extent was created out of the Bog Turtle dataset to limit the analysis to their geographic range and the extent was used to crop the variables. Finally, a raster layer was created with a prediction based on the model object. To view that information on a map you can use the gplot function which is a wrapper for ggplot from the rasterVis package.

gplot(presence_model) + 
  geom_polygon(data = us, aes(x= long, y = lat, group = group),
               fill = "gray", color="black") +
  geom_raster(aes(fill=value)) +
  geom_polygon(data = us, aes(x= long, y = lat, group = group),
               fill = NA, color="black") +
  geom_point(data = turt_df, aes(x = x, y = y), size = 2, color = "#ffff6d", alpha = 0.5) +
  scale_fill_gradientn(colours=c("#490092","#b66dff","#ffb6db"), "Probability") +
  coord_fixed(xlim = c(-84,-72), ylim = c(34.8,46)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Probability of Bog Turtle Occurrence") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + theme(legend.position = "right") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "lightblue"))

Now we have a beautiful map displaying the expected occurence of bog turtles according to the variables we tested. Because this analysis included “the kitchen sink” of variables, the model could be better specified by removing variables that are unrelated to the presence of Bog turtles. Additionally you can see areas where Bog Turtles has a high probability of occurrence but is not present. They may have occured there naturally and have now disappeared from habitat loss.

Let’s make this map interactive!

colors <- c("brown","yellow","darkgreen")

leaflet() %>% 
  addTiles() %>%
  addRasterImage(presence_model, colors = colors, opacity = 0.8) %>%
  addCircleMarkers(turt_df$x,
                   turt_df$y,
                   weight = 1,
                   color = "grey",
                   fillColor = "green",
                   fillOpacity = 0.7) %>%
  addMiniMap(position = 'topright',
             width = 100, 
             height = 100,
             toggleDisplay = FALSE) %>%
  addScaleBar(position = "bottomright")